In [ ]:
from os import path
import astropy.coordinates as coord
from astropy.io import fits
import astropy.units as u
import numpy as np
import matplotlib.pyplot as plt
from cycler import cycler
from itertools import cycle
plt.style.use('apw-notebook')
%matplotlib inline
import sqlalchemy
from sqlalchemy import func
from scipy.ndimage import gaussian_filter
from comoving_rv.db import Session, Base, db_connect
from comoving_rv.db.model import (Run, Observation, TGASSource, SimbadInfo, GroupToObservations,
SpectralLineInfo, SpectralLineMeasurement, RVMeasurement, Photometry)
In [ ]:
_color_cycler = cycler('color', ['#ea8b2c', '#7f5abf', '#1a9641', '#d7191c'])
color_cycler = cycle(_color_cycler)
In [ ]:
# base_path = '/Volumes/ProjectData/gaia-comoving-followup/'
base_path = '../../data/'
db_path = path.join(base_path, 'db.sqlite')
engine = db_connect(db_path)
session = Session()
In [ ]:
def get_abs_mag(mag, parallax, parallax_error):
# parallax in mas
SNR = parallax / parallax_error
dist = coord.Distance(1000. * (parallax/2 * (1 + np.sqrt(1 - 16/SNR**2)))**(-1) * u.pc)
mu = dist.distmod
M = mag - mu.value
return M
In [ ]:
tmass = fits.getdata('../../data/tgas_2mass_partial_j.fits.gz')
len(tmass)
In [ ]:
group_ids = session.query(Observation.group_id).join(Run).filter(Run.name == 'mdm-spring-2017')\
.filter((Observation.group_id != None) &
(Observation.group_id != 0) &
(Observation.group_id != 10))\
.distinct().all()
group_ids = [x[0] for x in group_ids]
len(group_ids)
In [ ]:
def get_color_mag(group_ids):
base_q = session.query(Observation).join(RVMeasurement).filter(RVMeasurement.rv != None)
color_mag = dict()
for gid in group_ids:
try:
gto = session.query(GroupToObservations).filter(GroupToObservations.group_id == gid).one()
obs1 = base_q.filter(Observation.id == gto.observation1_id).one()
obs2 = base_q.filter(Observation.id == gto.observation2_id).one()
except sqlalchemy.orm.exc.NoResultFound:
print('Skipping group {0}'.format(gid))
G_Js = []
M_Gs = []
for obs in [obs1, obs2]:
G = obs.tgas_source.phot_g_mean_mag
J = obs.tgas_source.J
if G is None or J is None:
break
G_Js.append(G-J)
M_Gs.append(get_abs_mag(G, obs.tgas_source.parallax,
obs.tgas_source.parallax_error))
else:
color_mag[gid] = {'G-J': np.array(G_Js), 'M_G': np.array(M_Gs)}
return color_mag
In [ ]:
color_mag_all = get_color_mag(group_ids)
In [ ]:
M_G_all = get_abs_mag(tmass['phot_g_mean_mag'], tmass['parallax'], tmass['parallax_error'])
G_J_all = tmass['phot_g_mean_mag'] - tmass['j_m']
xbins = np.arange(-0.1, 2.3+0.01, 0.02)
ybins = np.arange(-0.5, 8.5+0.01, 0.04)
# xbins = np.arange(-0.1, 2.3+0.01, 0.1)
# ybins = np.arange(-0.5, 8.5+0.01, 0.2)
H,xedges,yedges = np.histogram2d(G_J_all, M_G_all, bins=(xbins, ybins))
H = gaussian_filter(H, 1.5)
In [ ]:
derp = np.log(H.T+1.).ravel()
plt.hist(np.arctan(derp))
In [ ]:
plt.pcolormesh(xedges, yedges, np.arctan(np.log(H.T+1.)),
cmap='Blues', linewidth=0, rasterized=True,
alpha=1, edgecolors='None', vmin=0.2, vmax=2.75)
In [ ]:
def plot_cmd(color_mag, group_to_color=None, title='', markersize=4, ax=None):
if ax is None:
fig,ax = plt.subplots(1, 1, figsize=(6.5,7))
else:
fig = ax.figure
ax.pcolormesh(xedges, yedges, np.arctan(np.log(H.T+1.)),
cmap='Blues', linewidth=0, rasterized=True,
alpha=1, vmin=0.2, vmax=3.5)
if group_to_color is None:
group_to_color = dict()
for gid, d in color_mag.items():
color = group_to_color.get(gid, next(color_cycler)['color'])
l, = ax.plot(d['G-J'], d['M_G'], marker='', linewidth=1.,
linestyle='-', alpha=0.65, zorder=1, color=color)
group_to_color[gid] = color
ax.plot(d['G-J'], d['M_G'], marker='.',
linestyle='', alpha=1., color='k', zorder=10,
markersize=markersize)
ax.set_xlim(-0.1, 2.3)
ax.set_ylim(8.5, -0.5)
ax.set_xlabel('$G-J$ [mag]')
# ax.set_ylabel('$G - 5(\log\hat{d}-1)$ [mag]')
ax.set_ylabel('$M_G$ [mag]')
if title:
ax.set_title(title, fontsize=22)
return fig, group_to_color
In [ ]:
fig, group_to_color = plot_cmd(color_mag_all, title='All observed candidate comoving pairs')
fig.tight_layout()
# fig.savefig('sample_cmd.pdf', dpi=300)
In [ ]:
delta_M_G_all = []
for k in color_mag_all:
delta_M_G_all.append(abs(color_mag_all[k]['M_G'][0] - color_mag_all[k]['M_G'][1]))
In [ ]:
from astropy.table import Table
In [ ]:
tbl = Table.read('group_prob_dv.ecsv', format='ascii.ecsv')
In [ ]:
comoving = tbl['prob'] > 0.5
color_mag_genuine = get_color_mag(tbl['group_id'][comoving])
In [ ]:
delta_M_G_genuine = []
for k in color_mag_genuine:
delta_M_G_genuine.append(abs(color_mag_genuine[k]['M_G'][0] - color_mag_genuine[k]['M_G'][1]))
In [ ]:
fig,_ = plot_cmd(color_mag_genuine, group_to_color=group_to_color,
title='Probable comoving pairs')
fig.tight_layout()
fig.savefig('genuine_cmd.pdf', dpi=300)
In [ ]:
_,bins,_ = plt.hist(delta_M_G_all, bins='auto', normed=True, alpha=0.4)
plt.hist(delta_M_G_genuine, bins=bins, normed=True, alpha=0.4);
In [ ]:
from gwb.data import TGASData
In [ ]:
tgas = TGASData('../../../gaia-comoving-stars/data/stacked_tgas.fits')
In [ ]:
rave_star = Table.read('../../../gaia-comoving-stars/paper/t1-1-star.txt', format='ascii.csv')
In [ ]:
rave_tbl = Table.read('group_prob_dv_rave.ecsv', format='ascii.ecsv')
In [ ]:
rave_comoving = rave_tbl['prob'] > 0.5
color_mag = dict()
for gid in np.unique(rave_tbl['group_id'][rave_comoving]):
group = rave_star[rave_star['group_id'] == gid]
color_mag[gid] = {'G-J': [],
'M_G': []}
for row in group:
G = row['G']
J = row['J']
if G is None or J is None:
del color_mag[gid]
break
i1 = np.where(tgas._data['source_id'] == row['tgas_source_id'])[0][0]
star = tgas[i1]
M_G = get_abs_mag(G, star.parallax.value, star.parallax_error.value)
color_mag[gid]['G-J'].append(G - J)
color_mag[gid]['M_G'].append(M_G)
if gid in color_mag and len(color_mag[gid]['G-J']) < 2:
print('deleting')
del color_mag[gid]
for gid in np.unique(rave_tbl['group_id'][rave_comoving]):
if gid not in color_mag: continue
for k in color_mag[gid].keys():
color_mag[gid][k] = np.array(color_mag[gid][k])
In [ ]:
fig,_ = plot_cmd(color_mag, title='RAVE probable comoving pairs')
fig.savefig('genuine_cmd_rave.pdf', dpi=300)
Non-comoving pairs
In [ ]:
color_mag = dict()
for gid in tbl['group_id'][~comoving]:
group = session.query(Observation).join(Run).filter(Run.name == 'mdm-spring-2017')\
.filter(Observation.group_id == gid).all()
if len(group) > 2:
continue
color_mag[gid] = {'G-J': [],
'M_G': []}
for member in group:
src = member.tgas_source
G = src.phot_g_mean_mag
J = src.J
if G is None or J is None:
del color_mag[gid]
break
M_G = get_abs_mag(G, src.parallax, src.parallax_error)
color_mag[gid]['G-J'].append(G - J)
color_mag[gid]['M_G'].append(M_G)
if gid in color_mag and len(color_mag[gid]['G-J']) < 2:
del color_mag[gid]
for gid in tbl['group_id'][~comoving]:
if gid not in color_mag: continue
for k in color_mag[gid].keys():
color_mag[gid][k] = np.array(color_mag[gid][k])
In [ ]:
fig,ax = plt.subplots(1, 1, figsize=(6,6))
ax.pcolormesh(xedges, yedges, np.log(H.T+1.),
cmap='Blues', linewidth=0, rasterized=True)
for gid, d in color_mag.items():
if gid not in group_color_map:
continue
color = group_color_map[gid]
ax.plot(d['G-J'], d['M_G'], marker='', linewidth=1.,
linestyle='-', alpha=0.65, zorder=1, color=color)
ax.plot(d['G-J'], d['M_G'], marker='.',
linestyle='', alpha=1., color='k', zorder=10, markersize=3)
ax.set_xlim(-0.1, 2.3)
ax.set_ylim(8.5, -0.5)
ax.set_xlabel('$G-J$ [mag]')
# ax.set_ylabel('$G - 5(\log\hat{d}-1)$ [mag]')
ax.set_ylabel('$M_G$ [mag]')
ax.set_title('All observed comoving pairs')
fig.tight_layout()
# fig.savefig('genuine_cmd.pdf', dpi=300)
In [ ]:
interesting_group_ids = np.zeros(3, int)
In [ ]:
far_group_ids = [x[0] for x in session.query(Observation.group_id).join(TGASSource).filter(TGASSource.parallax < 10).all()]
In [ ]:
for gid, d in color_mag_genuine.items():
if (d['G-J']>1.1).any() and (d['G-J']<0.5).any() and np.all(d['M_G'] < 3):
print(gid)
interesting_group_ids[0] = gid
assert gid in far_group_ids
In [ ]:
obs1, obs2 = session.query(Observation).filter(Observation.group_id == 1500).all()
obs1.icrs().separation_3d(obs2.icrs()), obs1, obs2
In [ ]:
for gid, d in color_mag_genuine.items():
if ((d['G-J'][d['M_G'].argmin()] - d['G-J'][d['M_G'].argmax()]) > 0.2 and
np.all(d['G-J'] > 1) and
abs(d['M_G'][1]-d['M_G'][0]) > 1.5 and gid in far_group_ids):
print(gid)
interesting_group_ids[1] = gid
break
In [ ]:
obs1, obs2 = session.query(Observation).filter(Observation.group_id == interesting_group_ids[1]).all()
obs1.icrs().separation(obs2.icrs()), obs1, obs2
In [ ]:
(obs1.tgas_source.parallax, obs1.tgas_source.parallax_error), obs2.tgas_source.parallax
In [ ]:
delta_M_G_ = dict()
for gid, d in color_mag_genuine.items():
delta_M_G_[gid] = abs(d['M_G'][1]-d['M_G'][0])
max(list(delta_M_G_.items()), key=lambda x: x[1])
interesting_group_ids[2] = 1515
# assert interesting_group_ids[2] in far_group_ids
In [ ]:
obs1, obs2 = session.query(Observation).filter(Observation.group_id == 1515).all()
obs1.icrs().separation_3d(obs2.icrs()), obs1, obs2
In [ ]:
obs1.tgas_source.parallax, obs2.tgas_source.parallax
In [ ]:
interesting_group_ids
In [ ]:
interesting_color_mag = dict([(gid, color_mag_genuine[gid]) for gid in interesting_group_ids])
interesting_names = dict()
for gid in interesting_group_ids:
obs1, obs2 = session.query(Observation).filter(Observation.group_id == gid).all()
interesting_names[gid] = [obs1.simbad_info.preferred_name,
obs2.simbad_info.preferred_name]
In [ ]:
fig,_ = plot_cmd(interesting_color_mag, group_to_color=group_to_color,
title='Highlighted comoving pairs',
markersize=8)
ax = fig.axes[0]
offsets = {1515: [[0.1,0.25],
[-0.5,0.55]],
1500: [[-0.5,-0.55],
[-0.3,-0.2]],
1229: [[-.65,0.45],
[-0.75,0.6]]}
for gid in interesting_group_ids:
for i in range(2):
ax.text(interesting_color_mag[gid]['G-J'][i] + 0.02 + offsets[gid][i][0],
interesting_color_mag[gid]['M_G'][i] - 0.10 + offsets[gid][i][1],
interesting_names[gid][i], fontsize=20, zorder=10)
ax.arrow(x=1.1, y=2.3, dx=0.135, dy=0.37,
head_width=0.03, head_length=0.1,
linestyle='-', linewidth=1, color='#777777')
fig.savefig('highlighted_cmd.pdf', dpi=300)
In [ ]:
fig, axes = plt.subplots(1, 2, figsize=(11.25, 6.3), sharex=True, sharey=True)
fig,_ = plot_cmd(color_mag_genuine, group_to_color=group_to_color,
title='Probability $> 0.5$', ax=axes[0])
# highlighted
ax = axes[1]
fig,_ = plot_cmd(interesting_color_mag, group_to_color=group_to_color,
title='Highlighted comoving pairs',
markersize=8, ax=ax)
offsets = {1515: [[0.1,0.25],
[-0.5,0.55]],
1500: [[-0.5,-0.55],
[-0.3,-0.2]],
1229: [[-.65,0.45],
[-0.75,0.6]]}
for gid in interesting_group_ids:
for i in range(2):
ax.text(interesting_color_mag[gid]['G-J'][i] + 0.02 + offsets[gid][i][0],
interesting_color_mag[gid]['M_G'][i] - 0.10 + offsets[gid][i][1],
interesting_names[gid][i], fontsize=20, zorder=10)
ax.arrow(x=1.1, y=2.3, dx=0.135, dy=0.37,
head_width=0.03, head_length=0.1,
linestyle='-', linewidth=1, color='#777777')
ax.set_ylabel('')
fig.tight_layout()
fig.savefig('genuine_highlighted_cmd.pdf', dpi=300)
In [ ]: